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Abstract 

To  understand  antiwear  phenomena  in  motor  engines  at  the  atomic  level  and  provide  evidence  in 
selecting  future  ashless  wear  inhibitors,  we  studied  the  thermal  stability  of  the  self-assembled  monolayer 
(SAM)  model  for  dithiophosphate  (DTP)  and  dithiocarbamate  (DTC)  molecules  on  the  iron  oxide 
siuface  using  molecular  djuiamics.  The  interactions  for  DTP,  DTC  and  FeaOs  are  evaluated  based  on  a 
force  field  derived  from  fitting  to  ab  initio  quantum  chemical  calculations  of  dimethyl  DTP  (and  DTC) 
and  Fe(OH)2(HaO)3-DTP  (DTC)  clusters.  MD  simulations  at  constant-NPT  are  conducted  to  asses 
relative  thermal  stabilities  of  the  DTP  and  DTC  with  different  pendant  groups  (n-propyl,  i-propyl,  n- 
pentyl,  and  i-pentyl).  To  investigate  frictional  process,  we  employ  a  steady  state  MD  method,  in  which 
one  of  the  Fe203  slabs  maintmned  at  a  constant  linear  velocity.  We  obtain  the  time  averaged  normal 
and  frictional]  forces  from  the  interatomic  forces.  Then,  we  calculated  the  friction  coefficient  at  the 
interface  between  SAMs  of  DTP  and  the  confined  lubricant,  hexadecane,  to  assess  the  shear  stability  of 
DTPs  with  different  pendant  groups. 


1  Introduction 

Zinc  dithiophosphates  (ZnDTP)  have  been  used  for  several  decades'  to  reduce  wear  to  acceptable  levels  in 
modern  engines.  Although  DTP  is  the  most  effective  and  economic  antiwear  agent  used,  there  is  environ¬ 
mental  concern  due  to  the  level  of  phosphorous  content  in  the  formulation.  Finding  a  ZnDTP  replacement 
with  lower  phosphorus  content  and  providing  the  necessary  wear  protection  to  eiigine  surfaces  is  very  impor¬ 
tant  for  the  oil  industry.  One  might  consider  other  lubricant  additives  which  are  presently  used  to  control 
friction,  deposits,  corrosion,  oxidation,  and  rust  (for  a  review,  see  reference  2).  However,  the  number  of  pos¬ 
sible  chemical  compounds  is  enormous,  making  an  empirical  search  laborious  and  expensive.  It  should  also 
be  noted  that  a  single  complete  engine  test  for  a  new  wear  inhibitor  costs  around  $150,000.  Another  serious 
impediment  to  this  search  is  that  the  factors  underlying  good  wear  performance  are  not  well  understood. 

Using  atomistic  level  modelling  and  molecular  dynamics  (MD)  simulation  to  provide  guidance  in  pri¬ 
oritizing  new  materials  for  experimental  testing  is  a  viable  alternative.  For  instance,  our  earlier  work  led 
to  a  self-assembled  monolayer  (SAM)  model,  where  we  found  the  calculated  cohesive  energy  of  the  wear 
inhibitor  correlates  inversely  with  wear  performance.^''* 

The  purpose  of  a  lubricant  is  to  provide  a  protective  coating  to  the  solid  surfaces.  This,  in  turn  prevents 
formation  of  £idhesive  junctions,  reduces  frictional  energy  losses  by  acting  as  an  interfacial  layer  of  low 
shear  strength.  Our  understanding  of  the  friction  of  two  solid  surfaces  sliding  against  each  other  with  a 
thin  liquid  film  between  them  has  been  greatly  advanced  over  the  last  ten  years.  The  work  on  a  variety  of 
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different  surfaces,  such  as  mica,  silica,  metal  oxide,  and  surfactant  monolayer  surfaces,  sliding  across  different 
lubricants  have  shown  that  the  properties  of  fluids  confined  to  the  nanometer  scale  are  very  different  from 
the  bulk  properties.  For  example,  the  “effective”  viscosity  of  a  thin  lubricant  film  is  often  much  higher  than 
bulk  viscosity. 

In  this  paper,  extending  our  earlier  work  on  the  SAM  model  of  wear  inhibitors®’"*,  we  calculated  the 
adsorption  energies  and  coverage  for  DTP  and  DTCs  at  500  K.  We  found  that  for  DTP  the  SAM  model 
of  adsorption  energies  correlate  (inversely)  with  wear  obtained  from  engine  tests.  We  also  present  results 
of  molecular  dynamics  simulations  of  shear  stability  of  these  SAM  in  the  presence  of  lubricant  (normal- 
hexadecane).  The  lubricant  is  confined  by  two  iron-oxide  solid  surfaces  on  which  DTP  molecules  are 
chemisorbed  to  form  a  self-assembled  monolayer.  We  have  carried  out  the  simulations  for  three  different 
DTP  molecules:  isopropyl  (i-Pr),  isobutyl  (i-Bu),  and  phenyl  (Phe)  to  investigate  the  effect  of  the  SAMs 
on  the  dynamic  properties  of  the  lubricant. 

2  Simulation  Results  and  Discussion 

2.1  The  Force  Fields  for  DTP,  DTC  and  Iron  Oxide. 

We  modified  the  universal  force  field  (UFF)®  parameters  to  fit  the  ab  initio  calculations,  (geometry,  bond 
energy,  and  vibrational  frequency  from  HF  calculations  on  clusters  containing  DTP  bonded  to  the  metal) 
for  this  work. 

The  total  potential  energy  of  the  system  is  given  by 

E  =  E„ai  Efit  (1) 

The  valence  interactions,  Evai  =  Ebond  +  EangU  +  Etortion  +  Einvertion,  consist  of  bond  stretch  (Ebond), 
bond-angle  bend  (Eangie),  dihedral  angle  torsion  (Etorsion),  and  inversion  (Einvenion)  terms.  For  Ebond  we 
use  Morse  terms,  for  angle  terms  we  use  the  dihedral  cosine  harmonic  terms.  For  Etoraim  we  use  Etoraim  — 
2^(1+  cos  3<^)  and  for  the  inversion  energy,  we  have  used  the  harmonic  expression:  Einu  =  Ki  (ip  —  if>o)^- 

The  nonbonded  interactions,  Enb,  consist  of  vdw  (Evdw)  and  electrostatic  {Eq)  terms,  for  Evdw  we 
use  exponential-6  and  Morse  terms.  The  equilibrium  bond  lengths  (ilo)  and  angles  (l9o)  were  fitted  to 
the  ab  initio  geometries.  The  force  constants  (i.e.  Kb,  K$,  etc.)  were  fitted  to  the  ab  initio  vibrational 
frequencies  using  the  biased  Hessian  method.®  The  torsional  barriers  (Vm)  were  adjusted  to  reproduce  the 
conformational  energy  differences  from  the  HF  calculations.  For  metal-DTP  (-DTC)  systems  the  Fe-S 
interaction  was  treated  as  a  nonbonded  interaction  just  as  for  the  Zn-S  interactions  of  ZnDTP  systems.® 
That  is,  they  were  described  with  vdw  and  electrostatic  terms  only  (no  V2ilence  terms).  This  allows  DTP 
(or  DTC)  to  dissociate  smoothly  from  the  metal  surfaces.  The  further  details  and  parameters  are  given  in 
References  3  and  4. 

2.2  Self  Assembled  Monolayer  and  Adsorption  at  Different  Coverage 

The  simulation  results  presented  below  are  mainly  from  two  types  of  calculations:  one  from  minimized 
structures  and  the  other  from  MD  simulations  at  500  K  using  NVT  ensemble  molecular  dynamics.  The  run 
times  of  simulation  vary  between  10  ps  and  20  ps.  The  integeration  time  step  is  1  fs. 

To  determine  the  adsorption  energy,  we  first  calculate  the  cohesive  energy  as  a  function  of  the  number 
of  molecules  adsorbed, 

Ecoh  —  {Eaurf+SAM  Eturf  N  X  Ejnol^  (2) 

where  E,urf+SAM,  Eaurf,  and  Emoi  are  the  total  energy  of  the  Fe203  surface  plus  SAM,  that  of  the 
Fe203  surface  only,  and  that  of  the  SAM  molecule  in  the  gas  phase,  respectively.  The  adsorption  energy 
per  molecule,  Ead,  is  then  obtained  from  the  ration  of  the  cohesive  energy  to  the  number  of  adsorbate 
molecules. 

The  results  for  DTP  are  summarized  in  Figure  1,  (a)  from  the  minimized  structures,  and  (b)  from  the 
dynamics  simulations.  Clearly,  the  adsorption  energy  curves  in  Figure  1  are  quite  fla.t  when  the  number 
of  adsorbate  molecules  is  less  than  8  (i.e.,  <  1/2  coverage).  The  errors  in  adsorption  energy  under  our 
simulation  conditions  are  estimated  to  be  1  kcal/mol/molecule.  Within  this  error  the  adsorption  energies 
up  to  1/2  Coverage  can  be  considered  more  or  less  a  constant.  After  1/2  coverage,  the  adsorption  energies 
begin  to  drop.  The  considerable  reduction  in  adsorption  energy  indicates  that  the  ligand-effect  due  to 
interaction  between  neighboring  adsorbed  molecules  become  significant  when  the  surface  coverage  exceeds 
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Figure  1:  Adsorption  energy  as  a  function  of  surface  coverage  for  DTP  molecules:  (a)  from  minimized 
structures;  (b)  from  molecular  dynamics  simulations  at  500  K.  The  adsorption  energies  remain  almost  a 
constant  up  to  1/2  coverage  (i.e.,  8  adsorbate  molecules). 


Figure  2:  Correlation  of  SAM  adsorption  energy  with  wear  in  engine  wear  tests  for  commercial  DTP 
molecules.  The  error  in  the  adsorption  energy  is  1  kcal/mol/molecule.  The  straight  line  is  the  least-square 
fit  to  the  results  of  both  minimization  and  dynamics  simulation. 
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Figure  3:  Friction  simulation  setup;  2-D  periodic  boundary  conditions  imposed  in  x-  and  y-directions.  The 
top  layer  is  driven  at  a  constant  velocity  in  x-direction. 


1/2.  Similar  results  are  obtained  for  DTC.  Here,  the  decrease  in  adsorption  energy  for  these  DTC  molecules 
is  more  gradual,  and  it  does  not  start  exactly  at  1/2  coverage  in  some  cases.  For  example,  for  iso-heptyl  (iC?) 
and  normal-heptyl  (nCr),  the  adsorption  energies  begin  to  decrcEise  at  3/8  coverage.  However,  the  molecular 
dynamics  simulations  of  the  monolayers  of  these  DTC  molecules  show  an  enhancement  in  adsorption  energy 
as  the  critical  surface  coverage  increases  to  1/2. 

In  Figure  2,  we  plotted  the  adsorption  energy  at  1/2  coverage  against  the  measured  maximum  wear 
in  engine  wear  tests^  for  the  three  commercial  R  groups  of  DTP.  Both  the  minimization  and  dynamics 
simulation  results  are  shown.  The  straight  line  in  the  figure  is  the  least-square  fit  to  the  both  sets  of  data. 
We  see  that  the  adsorption  energies  correlate  with  the  wear  data.  A  stronger  adsorbed  monolayer  provides 
a  smaller  wear  and,  hence,  a  better  protection  to  the  surface.  This  correspondence  between  the  adsorption 
energy  and  wear  presents  a  basis  for  predicting  the  antiwear  performance  of  other  materials  prior  to  engine 
tests. 

2.3  Friction  and  Shear  Stability  of  SAM-Lubricant  Interface 

To  simulate  the  the  shear  stability  and  the  friction  between  the  SAM  of  wear  inhibitors  and  lubricant,  we 
considered  the  model  system  shown  in  Figure  3.  The  fluid  portion,  consisting  of  32  molecules  of  CieH^^ 
is  sandwiched  between  two  planar  walls.  Each  wall  is  made  of  an  hexagonal  iron-oxide  (0001)  surface  and 
a  monolayer  of  DTP  molecules;  the  binding  between  the  surface  and  SAM  is  through  chemisorption.  The 
geometrical  parameters  for  this  set  up  are,  D  =  2  nm,  H  =  4.4-4. 5  nm,  and  h=0.6-0.7  nm  depending  on 
the  type  of  SAM  used.  L,  is  around  8.5  nm,  the  actual  value  is  determined  by  the  energy  minimization  of 
the  entire  system  for  each  case. 

Simulations  aimed  at  generating  planar  Coueutte  flow  by  moving  the  top  wall  at  a  constant  velocity 
V  along  the  X  axis  while  holding  the  lower  wall  fixed.  In  the  results  presented  below,  v  =  1  A/ps.  2- 
dimensional  periodic  EC’s  are  imposed  along  the  X  and  Y  directions.  In  all  simulations,  the  time  step  is 
1  fs.  The  duration  of  runs  are  200  ps.  The  trajectories  are  stored  every  0.05  ps  for  analysis. 

In  order  to  capture  the  stick-slip  motion,  we  used  the  displacement  over  a  5  ps  interval  to  derive  the 
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Figure  4:  Velocity  profile  for  a)  i-Pr,  b)  i-Bu  and  c)  Phe.  The  crosses  are  the  average  streaming  velocity  for 
each  molecule,  the  solid  line  is  the  average  over  the  number  of  molecules 


streaming  velocity  of  the  molecules.  The  streaming  velocity  of  each  molecule  is  calculated  from  a  linear 
least-square  fit.  This  procedure  gives  a  smooth  streaming  velocity  profile  as  a  function  of  time  for  each 
molecule.  The  resulting  velocity  profiles  for  each  case  studied  (i-Pr,  i-Bu,  Phe)  are  shown  in  Figure  4.  The 
solid  line  in  each  panel  is  the  velocity  of  the  layer  (8-layers  constrcuted  along  Z-direction)  averaged  over 
time  and  over  the  number  of  the  molecules  in  that  layer.  The  crosses  on  each  panel  represent  the  streaming 
velocity  of  e2ich  molecule. 

It  is  apparent  that  the  different  characteristics  in  velocity  profile  have  resulted  from  the  difference  of  the 
strength  of  interactions  at  the  lubricant-SAM  interfaces.  Now  we  look  at  the  slip,  or  the  velocity  differences 
between  fluid  and  wall,  at  the  interfaces.  We  can  see  that  the  slip  for  i-Bu  is  larger  than  those  for  i-Pr  and 
Phe.  The  slip  at  the  wall-fluid  interfrice  is  one  of  issues  heavily  investigated  by  previous  simulations  because 
of  its  relevance  to  the  boundary  condition  in  fluid  dynamics.  A  broad  spectrum  of  boundary  conditions  was 
observed  ranging  from  pure  slip  to  complete  no-slip  (locking) .  The  degree  of  slip  at  the  boundary  depends 
on  several  interfacial  parameters  including  the  strength  of  the  liquid-solid  coupling,  the  thermal  roughness 
of  the  interface,  and  the  commensurability  of  wall  and  liquid  densities.  For  i-Pr,  the  profiles  indicate  a 
plug-flow  at  initial  times,  i.e.,  the  entire  fluid  is  dragged  along  by  the  top  wall.  This  over-shooting  is 
gradually  relaxed  in  later  times.  After  100  ps,  the  velocity  profile  approaches  a  steady-state.  This  kind  of 
flow  overshooting  is  much  less  profound  for  Phe,  and  it  is  inexistent  for  i-Bu. 

2.4  Properties:  Friction  Coefficient  and  Apparent  Viscosity 

The  normal  and  friction  forces  acting  on  the  fluid  due  to  the  interactions  from  the  top  and  bottom  walls  are 
calculated  for  each  run  (i-Pr,  i-Bu,  and  Phe)  to  determine  the  friction  coefficient  and  apparent  viscosity. 

It  seems  that  there  is  a  slight  increase  in  the  force  in  the  normal  direction  over  time.  This  may  be 
caused  by  the  increase  in  temperature  of  the  lubrication  during  the  dynamic.  In  these  simulations,  we  did 
not  implement  Langevin  bath  to  maintain  the  temperature  of  each  region.  So  that  we  could  also  calculate 
the  friction  force  from  the  slope  of  the  work  done  as  a  function  of  time.  The  two  methods  of  deriving  friction 
force,  from  direct  evaluation  or  from  the  slope  of  the  work,  give  the  consistent  results  on  the  friction  force. 
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(3) 


The  friction  coefficient  is  calculated  from  the  ratio  of  the  normal  force  and  friction  force: 


The  results  are  0.143  ±  0.001,  0.124  ±  0.00  and  0.135  ±0.001  nN  for  i-Pr,  i-Bu  and  Phe  cases,  respectively. 
It  is  seen  that  the  friction  coefficient  is  largest  for  i-Pr,  smallest  for  i-Bu,  and  falls  in  between  for  Phe.  But 
the  differences  may  be  insignificant  given  the  uncertainties. 

It  is  interesting  to  note  that  the  ordering  in  friction  coefficient  is  coincident  with  that  in  the  thickness 
of  SAM.  It  is  possible  that  a  thicker  SAM  diminishes  corrugation  of  the  underlying  surface  more  so  that 
the  friction  coefficient  is  reduced  when  the  DTP  molecules  with  longer  hydrocarbons  are  chemisorbed  to 
the  surface. 

In  experiments  the  effective  viscosity  is  calculated  from  the  me^ured  shear  stress  and  the  applied  shear 
rate  with  the  assumption  of  linear  velocity  profile  via 


a 

^  (dvldz) 


(4) 


where  er  is  the  shear  stress  and  {dvjdz)  is  the  shear  rate.  The  shear  stress  in  our  systems  are  obtained 
from  the  friction  force  divided  by  the  area  of  the  surface.  Since  the  velocity  profile  in  our  simulations  is  not 
linear  in  some  cases,  we  use  the  velocity  gradient  at  the  top  interface  to  estimate  the  effective  viscosity  of 
the  lubricant  close  to  the  top  interface.  For  a  direct  comparison,  the  friction  forces  of  different  SAMs  are 
normalized  to  the  same  normal  force.  The  apparent  viscosities  calculated  are,  2.68,  2.36,  2.41  mPa.s  (cP) 
for  i-Pr,  i-Bu  and  Phe  respectively. 

In  all  cases  we  studied,  the  surfaces  are  covered  by  a  chemisorbed  perfect  SAM  of  wear  inhibitors  and  a 
layer  of  lubricant.  In  our  simulations  no  wear  (large  friction  and  consequently  wear)  was  observed.  However, 
if  the  two  surfaces  were  in  direct  contact  in  some  region,  due  to  asperites,  the  high  friction  could  have  been 
strong  enough  to  tear  them  apart. 


3  Concluding  Remarks 

Assessment  of  a  new  material  demands  the  choice  of  a  suitable  criterion  that  can  distinguish  promising 
materials  from  those  with  little  potential.  For  wear  inhibitors,  the  key  parameters  that  we  considered  are  the 
adsorption  energy  at  minimization  and  at  elevated  temperatures,  and  the  shear  resistance  under  frictional 
conditions.  In  this  paper  we  presented  such  an  assessment  based  on  molecular  dynamics  simulations.  We 
observed  that  the  adsorption  energies  calculated  from  molecular  dynamics  simulations  at  500  K  correlate 
with  the  wear  data.  This  correlation  between  the  adsorption  energy  and  measured  wear  can  be  used  for 
predicting  the  antiwear  performance  of  other  materials  prior  to  engine  tests. 
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